####################################################################################
####################################################################################
# Appendix M: Tables M.57-M.60
####################################################################################
library(dplyr)
library(tidyr)
library(Hmisc)
library(lmtest)
library(multiwayvcov)
library(knitr)
library(stargazer)
library(ggplot2)
library(here)

dat = readRDS(here("data/dat_unrestricted.rds"))
datb = readRDS(here("data/dat_restricted.rds"))

####################################################################################
## Unrestricted

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Bushenyi","Ibanda","Jinja","Mbale","Mpigi","Pallisa","Mbarara","Shema",
         "village_distance_HF","village_distance_HOSP","village_distance_road","village_distance_electricity", "share2006")
dat[, paste0("c_", cols)] = lapply(dat[, cols], function(x) (x - mean(x))  )
covariates = c("treatment*c_village_distance_HF + treatment*c_village_distance_HOSP + treatment*c_village_distance_road + treatment*c_village_distance_electricity + 
                treatment*c_sex + c_age*treatment + c_education_level*treatment + treatment*c_share2006 + 
                c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment",
                "c_Arua*treatment + c_Bushenyi*treatment + c_Ibanda*treatment + c_Jinja*treatment + c_Mbale*treatment + c_Mpigi*treatment + c_Pallisa*treatment + c_Mbarara*treatment")

####################################################################################
## Restricted

# Implement SWE of ATE: Demean all covariates, fully interact
cols = c("sex","age","education_level","Arua","Ibanda","Pallisa","Shema",
         "village_distance_HF","village_distance_HOSP","village_distance_road","village_distance_electricity", "share2006")
datb[, paste0("c_", cols)] = lapply(datb[, cols], function(x) (x - mean(x))  )
covariatesb = c("treatment*c_village_distance_HF + treatment*c_village_distance_HOSP + treatment*c_village_distance_road + treatment*c_village_distance_electricity + 
                treatment*c_sex + c_age*treatment + c_education_level*treatment + treatment*c_share2006 + 
                c_Arua*treatment + c_Ibanda*treatment + c_Pallisa*treatment", 
               "c_Arua*treatment + c_Ibanda*treatment + c_Pallisa*treatment")


####################################################################################
####################################################################################
## Political Credit ####
# credit_i_index = Influence index
# perform_i_past_i_index = Past performance index
# perform_i_index = Current performance index
# ca_i_index = Full index index

## H3: Full government credit index

# Unrestricted
for (i in names(dat[,c("ca_gov_index","ca_pres_index","ca_gad_index","ca_gak_index","ca_mp_index","ca_lc_index","ca_dc_index","ca_ngo_index")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

# Restricted
for (i in names(datb[,c("ca_gov_index","ca_pres_index","ca_gad_index","ca_gak_index","ca_mp_index","ca_lc_index","ca_dc_index","ca_ngo_index")])){ # Dependent Variables
  for (j in covariatesb) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("bm",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = datb))
    # Output clustered SEs (county)
    assign(x = paste("bc",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = datb),
                            cluster.vcov(lm(as.formula(model), data = datb), datb$village_id)))
  }
}

# stargazer(m.ca_lc_index.t, m.ca_dc_index.t, m.ca_mp_index.t, m.ca_pres_index.t, m.ca_gad_index.t, m.ca_gak_index.t,
#           
#           se = list(c.ca_lc_index.t[,2], c.ca_dc_index.t[,2], c.ca_mp_index.t[,2], c.ca_pres_index.t[,2], 
#                     c.ca_gad_index.t[,2], c.ca_gak_index.t[,2]),
#           
#           p = list(c.ca_lc_index.t[,4], c.ca_dc_index.t[,4], c.ca_mp_index.t[,4], c.ca_pres_index.t[,4], 
#                    c.ca_gad_index.t[,4], c.ca_gak_index.t[,4]),
#           
#           keep = c("treatment$"), 
#           
#           order = c("$treatment$"), covariate.labels=c("Treatment"),
#           
#           type = "latex", out = "tables/credit_1_urc.tex",
#           label = "tab:credit_1_urc", column.sep.width = "1pt", table.placement = "!ht", dep.var.caption = "",
#           keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
#           title = "Effect of LG CHP on Credit for Political Actors",
#           notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", notes.align = "l", notes.append = F, notes.label = "",
#           column.labels = c("Local", "District","MP", "President","District","National\\\\ & Councilors & Chair & & & Agency & Agency"))

####################################################################################
####################################################################################
## Table M.57

stargazer(m.ca_lc_index.c, m.ca_dc_index.c, m.ca_mp_index.c, m.ca_pres_index.c, m.ca_gad_index.c, m.ca_gak_index.c, 
          
          se = list(c.ca_lc_index.c[,2], c.ca_dc_index.c[,2], c.ca_mp_index.c[,2], c.ca_pres_index.c[,2], 
                    c.ca_gad_index.c[,2], c.ca_gak_index.c[,2]),
          
          
          p = list(c.ca_lc_index.c[,4], c.ca_dc_index.c[,4], c.ca_mp_index.c[,4], c.ca_pres_index.c[,4], 
                   c.ca_gad_index.c[,4], c.ca_gak_index.c[,4]),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/credit_1_urnc.tex",
          label = "tab:credit_1_urnc", column.sep.width = "1pt", table.placement = "!ht", dep.var.caption = "",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Credit for Political Actors (No covariates)",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Local", "District","MP", "President","District","National\\\\ & Councilors & Chair & & & Agency & Agency"))


# stargazer(bm.ca_lc_index.t, bm.ca_dc_index.t, bm.ca_mp_index.t, bm.ca_pres_index.t, bm.ca_gad_index.t, bm.ca_gak_index.t, 
#           
#           se = list(bc.ca_lc_index.t[,2], bc.ca_dc_index.t[,2], bc.ca_mp_index.t[,2], bc.ca_pres_index.t[,2], 
#                     bc.ca_gad_index.t[,2], bc.ca_gak_index.t[,2]),
#           
#           
#           p = list(bc.ca_lc_index.t[,4], bc.ca_dc_index.t[,4], bc.ca_mp_index.t[,4], bc.ca_pres_index.t[,4], 
#                    bc.ca_gad_index.t[,4], bc.ca_gak_index.t[,4]),
#           
#           keep = c("treatment$"), 
#           
#           order = c("$treatment$"), covariate.labels=c("Treatment"),
#           
#           type = "latex", out = "tables/credit_1_rc.tex",
#           label = "tab:credit_1_rc", column.sep.width = "1pt", table.placement = "!ht", dep.var.caption = "",
#           keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
#           title = "Effect of LG CHP on Credit for Political Actors (Restricted)",
#           notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", notes.align = "l", notes.append = F, notes.label = "",
#           column.labels = c("Local", "District","MP", "President","District","National\\\\ & Councilors & Chair & & & Agency & Agency"))

####################################################################################
####################################################################################
## Table M.58

stargazer(bm.ca_lc_index.c, bm.ca_dc_index.c, bm.ca_mp_index.c, bm.ca_pres_index.c, bm.ca_gad_index.c, bm.ca_gak_index.c,
          
          se = list(bc.ca_lc_index.c[,2], bc.ca_dc_index.c[,2], bc.ca_mp_index.c[,2], bc.ca_pres_index.c[,2],
                    bc.ca_gad_index.c[,2], bc.ca_gak_index.c[,2]),
          
          
          p = list(bc.ca_lc_index.c[,4], bc.ca_dc_index.c[,4], bc.ca_mp_index.c[,4], bc.ca_pres_index.c[,4],
                   bc.ca_gad_index.c[,4], bc.ca_gak_index.c[,4]),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/credit_1_rnc.tex",
          label = "tab:credit_1_rnc", column.sep.width = "1pt", table.placement = "!ht", dep.var.caption = "",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Credit for Political Actors (Restricted; No covariates)",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Local", "District","MP", "President","District","National\\\\ & Councilors & Chair & & & Agency & Agency"))

rm(list =  grep("^(c|m)\\.", ls(), value = T))
rm(list = grep("^(bc|bm)\\.", ls(), value = T))

## H3: President Components

# Unrestricted
for (i in names(dat[,c("ca_pres_index","performance_health_pres_st","performance_pres_st", "power_pres_st")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

# Restricted
for (i in names(datb[,c("ca_pres_index","performance_health_pres_st","performance_pres_st", "power_pres_st")])){ # Dependent Variables
  for (j in covariatesb) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("bm",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = datb))
    # Output clustered SEs (county)
    assign(x = paste("bc",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = datb),
                            cluster.vcov(lm(as.formula(model), data = datb), datb$village_id)))
  }
}

# stargazer(m.ca_pres_index.t, bm.ca_pres_index.t, m.power_pres_st.t, bm.power_pres_st.t, m.performance_health_pres_st.t, 
#           bm.performance_health_pres_st.t,m.performance_pres_st.t, bm.performance_pres_st.t, 
#           
#           se = list(c.ca_pres_index.t[,2], bc.ca_pres_index.t[,2], c.power_pres_st.t[,2], bc.power_pres_st.t[,2], 
#                     c.performance_health_pres_st.t[,2], bc.performance_health_pres_st.t[,2],c.performance_pres_st.t[,2], bc.performance_pres_st.t[,2] ),
#           
#           p = list(c.ca_pres_index.t[,4], bc.ca_pres_index.t[,4], c.power_pres_st.t[,4], bc.power_pres_st.t[,4], 
#                    c.performance_health_pres_st.t[,4], bc.performance_health_pres_st.t[,4],c.performance_pres_st.t[,4], bc.performance_pres_st.t[,4] ),
#           
#           keep = c("treatment$"), 
#           
#           order = c("$treatment$"), covariate.labels=c("Treatment"),
#           
#           type = "latex", out = "tables/credit_pres_c.tex",
#           label = "tab:credit_pres_c", column.sep.width = "1pt", table.placement = "!ht",
#           keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
#           title = "Effect of LG CHP on Credit to the President",
#           notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
#           column.labels = c("Index","Power","Health","General"), column.separate = c(2,2,2,2),
#           add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

####################################################################################
####################################################################################
## Table M.59

stargazer(m.ca_pres_index.c, bm.ca_pres_index.c, m.power_pres_st.c, bm.power_pres_st.c, m.performance_health_pres_st.c, 
          bm.performance_health_pres_st.c,m.performance_pres_st.c, bm.performance_pres_st.c, 
          
          se = list(c.ca_pres_index.c[,2], bc.ca_pres_index.c[,2], c.power_pres_st.c[,2], bc.power_pres_st.c[,2], 
                    c.performance_health_pres_st.c[,2], bc.performance_health_pres_st.c[,2],c.performance_pres_st.c[,2], bc.performance_pres_st.c[,2] ),
          
          p = list(c.ca_pres_index.c[,4], bc.ca_pres_index.c[,4], c.power_pres_st.c[,4], bc.power_pres_st.c[,4], 
                   c.performance_health_pres_st.c[,4], bc.performance_health_pres_st.c[,4],c.performance_pres_st.c[,4], bc.performance_pres_st.c[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/credit_pres_nc.tex",
          label = "tab:credit_pres_nc", column.sep.width = "1pt", table.placement = "!ht",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Credit to the President (No Covariates)",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("Index","Power","Health","General"), column.separate = c(2,2,2,2),
          add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes", "No", "Yes")))

rm(list =  grep("^(c|m)\\.", ls(), value = T))
rm(list = grep("^(bc|bm)\\.", ls(), value = T))


####################################################################################
## Substitution ####
# CHP Contact
# 4 "Within the past one month" 3 "Within the past three months" 2 "Within the past six months" 1 "Within the past year" -1 "More than 1 year ago" 0 "Never" -99 "Don't Know"

# Unrestricted
for (i in names(dat[,c("satisfied_st","total_visits_st")])){ # Dependent Variables
  for (j in covariates) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("m",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = dat))
    # Output clustered SEs (county)
    assign(x = paste("c",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = dat),
                            cluster.vcov(lm(as.formula(model), data = dat), dat$village_id)))
  }
}

# Restricted
for (i in names(datb[,c("satisfied_st","total_visits_st")])){ # Dependent Variables
  for (j in covariatesb) {
    model = paste(i,"~","treatment","+", j)
    
    # Run each model
    assign(x = paste("bm",i,substr(j,1,1), sep = "."), 
           value = lm(as.formula(model), data = datb))
    # Output clustered SEs (county)
    assign(x = paste("bc",i,substr(j,1,1),sep = "."), 
           value = coeftest(lm(as.formula(model), data = datb),
                            cluster.vcov(lm(as.formula(model), data = datb), datb$village_id)))
  }
}

# stargazer(m.satisfied_st.t, bm.satisfied_st.t, m.total_visits_st.t, bm.total_visits_st.t, 
#           
#           se = list(c.satisfied_st.t[,2], bc.satisfied_st.t[,2], c.total_visits_st.t[,2], 
#                     bc.total_visits_st.t[,2] ),
#           
#           p = list(c.satisfied_st.t[,4], bc.satisfied_st.t[,4], c.total_visits_st.t[,4], 
#                    bc.total_visits_st.t[,4] ),
#           
#           keep = c("treatment$"), 
#           
#           order = c("$treatment$"), covariate.labels=c("Treatment"),
#           
#           type = "latex", out = here("tables/subset_controls2/substitute_1_c.tex"),
#           label = "tab:substitute_1_c",column.sep.width = "1pt", table.placement = "!ht",
#           keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
#           title = "Effect of LG CHP on Perceptions and Use of VHTs",
#           notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
#           column.labels = c("VHT Satisfaction","VHT Use"), column.separate = c(2,2),
#           add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes")))

####################################################################################
####################################################################################
## Table M.60

stargazer(m.satisfied_st.c, bm.satisfied_st.c, m.total_visits_st.c, bm.total_visits_st.c, 
          
          se = list(c.satisfied_st.c[,2], bc.satisfied_st.c[,2], c.total_visits_st.c[,2], 
                    bc.total_visits_st.c[,2] ),
          
          p = list(c.satisfied_st.c[,4], bc.satisfied_st.c[,4], c.total_visits_st.c[,4], 
                   bc.total_visits_st.c[,4] ),
          
          keep = c("treatment$"), 
          
          order = c("$treatment$"), covariate.labels=c("Treatment"),
          
          type = "latex", out = "tables/substitute_1_nc.tex",
          label = "tab:substitute_1_nc",column.sep.width = "1pt", table.placement = "!ht",
          keep.stat = c("n"), dep.var.labels.include = F, no.space = T, model.numbers = T,
          title = "Effect of LG CHP on Perceptions and Use of VHTs (No covariates)",
          notes = "Standard errors are clustered at the village level. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01", dep.var.caption = "", notes.align = "l", notes.append = F, notes.label = "",
          column.labels = c("VHT Satisfaction","VHT Use"), column.separate = c(2,2),
          add.lines = list(c("Restricted", "No", "Yes", "No", "Yes", "No", "Yes")))

rm(list =  grep("^(c|m)\\.", ls(), value = T))
rm(list = grep("^(bc|bm)\\.", ls(), value = T))
